classdef AlkaliMetal < Atom.Atom
    %ALKALIMETAL Summary of this class goes here
    %   Detailed explanation goes here
    properties
        parameters         % alkali metal parameters
        hasSpin = 1

        I                  % alkali metal nuclear spin
        S = 0.5

        a                 % a:the atomic spin quantum number correspongding to a=F=I-J
        b                  % b:the atomic spin quantum number correspongding to a=F=I+J

        J=[0.5, 0.5, 1.5] % total angular momentum  of electron-S=1/2,J=1/2,J=3/2

        gI                % dimension of the nuclear spin space
        gJ                % dimension of the electron spin space with total angular momentum J

        dim               % dimension of the direction product space {|J,M_J>|I,M_I>}

        mat               % revalent matrix in the uncoupled basis
        matEigen
        eigen

        operator         % operators in Liouville space (used in master equation)

%        zeeman_terms

%         rho
    end

    methods
        function obj = AlkaliMetal(name, coil)
            obj@Atom.Atom(name);
            obj.parameters = Atom.AtomParameters(name);

            obj.I = obj.parameters.spin_I;
            obj.a = obj.I-0.5;
            obj.b = obj.I + 0.5;

            obj.gI = 2*obj.I+1;
            obj.gJ = 2*obj.J+1;
            obj.dim = obj.gI * obj.gJ;          % dimension of the space spanned by {|J,M_J>|I,M_I}

            obj.mat = obj.spinMatrix();

%             obj.rho = Algorithm.DensityMatrix(obj);
            if nargin > 1
                obj.set_eigen(coil);
            end

        end

        function obj = set_coil(obj, coil)
            obj.set_eigen(coil);
        end

        function obj = set_eigen(obj, coil)  % input (coil) is a 1X3 cell array corresponding to B_x,B_y,B_z.but only B_z is used.
            obj.eigen=Algorithm.Eigen(obj, coil{3});

            obj.matEigen.Imat = obj.eigen.transform(obj.mat.Imat);
            obj.matEigen.Smat = obj.eigen.transform(obj.mat.Smat);
            obj.matEigen.Fmat = obj.eigen.transform(obj.mat.Fmat);
            obj.matEigen.IS = obj.eigen.transform(obj.mat.IS);
            obj.matEigen.F2 = obj.eigen.transform(obj.mat.F2);
            obj.matEigen.mu = obj.eigen.transform(obj.mat.mu);
            for k=1:3
                obj.matEigen.zeeman_terms(:,:,k) = - obj.matEigen.mu{1}(:,:,k)*coil{k}.magB/(2*pi*h_bar)*1e-6; % energy in MHz
            end

%             obj.dipole();
            obj.operator = obj.spinOperator();
        end

        function res = energy_spectrum(obj, state, coil)
            coilx=coil{1};
            coily=coil{2};
            coilz=coil{3};
            res.magB = coilz.iter.dataList;
            res.energy = zeros(obj.dim(state), ...
                               coilz.iter.length);

            coilz.restart0();
            while coilz.iter.hasNext
                obj.set_eigen( {coilx, coily, coilz.move_forward} );
                res.energy(:, coilz.iter.cursor) ...
                    = obj.eigen.getEnergy(state);
            end
        end

        function S = mean_spin(obj, rho)
            s_mat = obj.matEigen.Smat{Atom.Subspace.GS};
            S = rho.mean(s_mat);
        end

        function v = mean_vector( obj, rho, space )
            if nargin < 3
                space = Atom.Subspace.GS;
            end

            i_mat = obj.matEigen.Imat{space};
            s_mat = obj.matEigen.Smat{space};
            f_mat = obj.matEigen.Fmat{space};
            v.I = rho.mean(i_mat);
            v.S = rho.mean(s_mat);
            v.F = rho.mean(f_mat);
        end

        function p = proj_op(obj, k, space)
            if nargin < 3
                space = Atom.Subspace.GS;
            end
            p = zeros( obj.dim(space) );
            p(k,k) = 1;
        end

    end

end
